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Coevolving and competing species or game-theoretic strategies exhibit rich and complex dynamics for which 
a general theoretical framework based on finite populations is still lacking. Recently, an explicit mean-field 
description in the form of a Fokker-Planck equation was derived for frequency-dependent selection with two 
strategies in finite populations based on microscopic processes [A. Traulsen, J. C. Claussen, and C. Hauert, 
Phys. Rev. Lett. 95, 238701 (2005)]. Here we generalize this approach in a twofold way: First, we extend the 
framework to an arbitrary number of strategies and second, we allow for mutations in the evolutionary process. 
The deterministic limit of infinite population size of the frequency dependent Moran process yields the adjusted 
replicator-mutator equation, which describes the combined effect of selection and mutation. For finite popula- 
tions, we provide an extension taking random drift into account. In the limit of neutral selection, i.e. whenever 
the process is determined by random drift and mutations, the stationary strategy distribution is derived. This 
distribution forms the background for the coevolutionary process. In particular, a critical mutation rate Uc is 
obtained separating two scenarios: above Uc the population predominantly consists of a mixture of strategies 
whereas below Uc the population tends to be in homogenous states. For one of the fundamental problems in 
evolutionary biology, the evolution of cooperation under Darwinian selection, we demonstrate that the analytical 
framework provides excellent approximations to individual based simulations even for rather small population 
sizes. This approach complements simulation results and provides a deeper, systematic understanding of coevo- 
lutionary dynamics. 

PACS numbers: 87.23.-n, 89.65.-s 87.23.Kg 



I. INTRODUCTION 

Darwinian evolution incorporates three basic processes: 
mutation, selection and random drift. In evolving popula- 
tions mutations create variation, which produce fitness differ- 
ences for selection to act upon and random drift accounts for 
stochastic effects arising due to the finite size of the popula- 
tion. The consideration of finitepopulations has a long tradi- 
tion in population genetics Ql |2 yl 0] with frequency inde- 
pendent selection. Often, finite size effects can be viewed as 
corrections to a continuum theory based on an infinite popula- 
tion mi^. Here, we analyze such finite population effects for 
frequency dependent selection. Under frequency dependent 
selection, or co-evolutionary dynamics, the fitness of each 
individual is affected by interactions with other members of 
the population. The corresponding mathematical framework 
is evolutionary game theory 0]. In the absence of muta- 
tions and in the deterministic limit of large well-mixed pop- 
ulations where individuals interact randomly, the replicator 
equation l^ fioll describes the change in frequency of the dif- 
ferent strategic types in the population. Including mutations 
in this limit leads to the replicator-mutator equation 
Interestingly, results can be quite different when turning to 
coevolutionary processes in finite populations based on the 
Moran process The frequency dependent Moran 

process is a stochastic birth-death process where an individ- 
ual is randomly selected for reproduction with a probability 
proportional to its fitness. The clonal offspring then replaces 
a randomly chosen individual in the population such that the 
population size remains constant. If there are only two types 
of individuals in the population, the evolutionary process can 
be mapped onto a random walk in one dimension IT^ . For 



more than two types the dynamics becomes significantly more 
difficult O- Further complications arise when individuals no 
longer interact randomly but, instead, spatial structure leads to 
limited local interactions. Traditionally, this is modeled by ar- 
ranging individuals on a spatial lattice or on more 
general networks [20^ 21^ 22^ 22j 2Ai 25 | and leads to new 
phenomena, which includes critical phase transitions j23l . and 
can affect the evolutionary process even in frequency indepen- 
dent settings |27, 28]. 

For two strategic types in well-mixed populations we estab- 
lished an explicit connection between microscopic stochastic 
processes in finite populations and the deterministic limit of 
infinite populations 12^ . The finite size effects are captured 
in the drift and diffusion terms of a Fokker-Planck equation 
where the diffusion term vanishes with 1/iV for increasing 
population sizes. In the present paper the Fokker-Planck equa- 
tion is generalized to an arbitrary number of strategies as well 
as to cover mutational changes of the strategic type. 

In Secinj the frequency dependent Moran process is gen- 
eralized to d strategies including mutations. The derivation of 
the corresponding Fokker-Planck equation is recalled in Sec. 
|ni]and the stationary distribution is derived in Sec. II VI This 
covers frequency independent and frequency dependent selec- 
tion. As applications, we consider the Prisoner's Dilemma and 
the Snowdrift game. The relation to dynamics in infinite pop- 
ulations is derived in Sec. |V| An example for an alternative 
microscopic process is provided in Sec. lVIl 

II. GENERALIZED MORAN PROCESS 

The frequency dependent Moran process has been intro- 
duced as a stochastic process among two types fisll . Although 
the extension to more strategies is straightforward fl^ . the 
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dynamics becomes significantly more difficult. In complete 
analogy to the Moran process among two individuals, an in- 
dividual is chosen at random proportional to fitness in every 
time step. This individual produces offspring which replaces 
a randomly chosen individual. In addition to the original fre- 
quency dependent Moran process, we include mutations, i.e. 
an individual of type k produces offspring of type j with prob- 
abiliy qjk, where '^'^=1 Ijk = 1 and d is the number of differ- 
ent types in the population. The special case of vanishing mu- 
tations is recovered by setting qjk = Sjk- Note that there are 
no restrictions on the frequency of mutations or on the states in 
which mutations occur For d > 2, the state space of the sys- 
tem is not longer a simple one-dimensional chain, but the dis- 
cretized simplex = {(ii, . . . ,id) & : *i = 

where ij is the number of individuals of type j and N is the 
total size of the population. The fitness ttj of an individual 
of type j is a linear combination of the payoff from interac- 
tions given by the entries of the d x d payoff matrix M with 
elements rrijk and a baseline fitness, which is set to one for 
convenience. Hence, 
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where self interactions are excluded by subtracting the payoff 
nijj. The selection intensity w determines the relative con- 
tributions of the baseline fitness and the interactions to the 
total fitness of an individual. For w ^ selection is weak 
and payoffs from the game represent small perturbations to 
the constant background fitness. For strong selection, w ^ 1, 
the influence of the background fitness vanishes and fitness is 
determined entirely by interactions. Since individuals are se- 
lected proportional to fitness, the fitness has to be a positive 
number. This can be guaranteed by payoff matrices with pos- 
itive entries only, rrijk > for all j and k, or by a sufficiently 
small w. 

Based on this fitness definition, the probability Tkj that an 
individual of type k is replaced by an individual of type j can 
be calculated. Tkj is the product of the probability that the 
individual chosen at random for elimination is of type k and 
the probability that an offspring of type j is produced. There 
are two possibilities to produce such an offspring: An indi- 
vidual of type j is chosen for reproduction with probability 
(Ncj)), where 4> = X]l=i T^mim/N, is the average pay- 
off in the population. It produces identical offspring (with 
probability qjj). The second possibility is that an individual 
of type I produces an offspring of type j due to mutations 
(with probability qij). Hence, the probability Tkj is given by 
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These transition probabilities describe the effect of selection 
and mutation in finite populations. The probabilistic nature of 
this coevolutionary process accounts for random drift. 

Since the state space of this system is no longer a one- 
dimensional chain for d > 2, many standard methods can no 
longer be applied. However, as in the case d ~ 2 a Fokker- 



Planck equation for the dynamics of the system can be derived 
for large populations. 

Although we are motivated by the generalized frequency 
dependent Moran process, the calculations following in Sees. 
Uni and IIVI are valid in a more general sense. In particular, 
they apply to any coevolutionary birth death process, such as 
the local update rule discussed in ll29ll . see Sec. lVIl 



in. FOKKER-PLANCK EQUATION 



In analogy to Ref. (291 for d = 2, a general Fokker-Planck 
equation for the probability density of strategies p{x) in an 
evolutionary game among d types can be derived. The con- 
stant population size of the generalized Moran process leads 
to a normalization of x, •'^j ~ ^- Thus, we have — 1 

independent variables. As shown in the Appendix, a Kramers- 
Moyal expansion yields the Fokker-Planck equation: 
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-p{x)bjk{x), 
(3) 

where x = {xi, . . . , Xd) with Xj = ij/N. The drift coeffi- 
cients are given by 



ak{x)=Y,Tjk{x)-Tkj{x) 



(4) 



and can be interpreted as an effective flow into state k: The 
first term is the probability that the number of k individuals 
increases, whereas the second term describes the probability 
for transitions from k to other states. For the diffusion matrix, 
we find 



bikix)= — 



-T,k{x)^ Tk, [x) + 5,k X T3I [x) + Tij {x) 



1=1 



(5) 

For two types, d = 2, and without mutations, the Fokker- 
Planck equation from l29il is recovered. Eq. (|3} describes 
the dynamics on the basis of changes in the probability dis- 
tribution. Traditionally, coevolutionary systems are often de- 
scribed by considering dynamical equations for the state of 
the system. General conclusions are then based on averages 
made over several realizations of the process. Here, the corre- 
spondence between both descriptions becomes clear, since the 
Fokker-Planck equation (|3ji corresponds to a stochastic dif- 
ferential equation |30, 31, 32]. The noise arises only from 
stochastic updating and is therefore not correlated in time. 
Hence, the Ito calculus can be applied to derive a Langevin 
equation describing the development of x. Equation (|3} cor- 
responds to the stochastic replicator-mutator equation 



Xk ^ akix) + y_] Ckj ix)Cj (<) 



(6) 



where Ckj is defined by J2'i=i Cki{x)cij{x) = bkj{x). Each 
element of the vector ^ is Gaussian white noise with unit vari- 
ance. The different elements are uncorrected, {^k{t)^j{s)) = 



3 



5kj5{t — s). Equation (|6} allows to approximate fluctuations 
arising from finite populations by Langevin terms that appear 
in a replicator equation. For any given payoff matrix, selec- 
tion pressure, population size and mutation rate, it provides 
a quantitative description of the fluctuations introduced by 
stochastic microscopic update processes. 



The drift and diffusion terms in the Fokker-Planck equation 
become 



a{x) = x{\ — x) 



ita{,x) ~ ttb{x) 



(1 - x)ttb{x) - xiia{.x) 



(13) 



IV. STATIONARY DISTRIBUTION 

The stationary distribution of strategies p*{x) can be de- 
rived from the Fokker-Planck equation which can be writ- 
ten in the form 



where the d 



p{x) = -VJ{x), (7) 
1 elements of the probabiUty current J{x) are 
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p{x)ak{x) - T^Yl -Q^.P'y^)^3k{x). 



(8) 



In general, the stationary solution p{x) = is equivalent to a 
probability current without sources. A special case of this is a 
probability current vanishing everywhere, i.e. Jk{x) ~ for 
all X. This leads to 
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If we exclude the degenerate cases with deihjk{x) = 0, the 
matrix bjk{x) can be inverted and Eq. (|9} reduces to 
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(10) 

A solution of this equation only exists if r(a;) is a gradient 
I30II32I1 . In that case, the stationary solution of Eq. (|3} can 
be obtained from a line integral between xo and x (which is 
independent of the path) 



p*{x) = A/'cxp 



/ r(y).dy 



(11) 



where M is a. normalization constant. However, in co- 
evolutionary systems the requirement that T{x) is a gradient 
is often not fulfilled. In that case, a stationary solution may 
still exist in which the ansatz of a vanishing probability cur- 
rent, J{x) = 0, is not valid although the divergence of J{x) 
remains zero, WJ{x) — 0. In these cases the stationary dis- 
tribution has to be derived by different means. 

A particularly simple case where these problems do not 
arise is d = 2, where only two types A and B are present. 
The state of the system can be described by the fraction x of 
A individuals. The remaining fraction 1 — x consists of B in- 
dividuals. Assuming a mutation rate u from ^ to i? as well as 
from B to A, the transition probabilities read 

XTTAix){l — m) + (1 — x)TrBix)u 



Tba{x) 
Tab [x) 
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xita{x)u + (1 — x)i:b{x){1 — u) 
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u xt:a{x) + (a; — 1)1: b{x) 
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(2a;- 1). 



The stationary solution of Eq. (|3} can now be computed as 
outlined above and is given by 



p*{x) = A/'exp 



(15) 



where V{y) = b-\y){2a{y) - b'{y)), b'iy) = j-^b{y) and 

M = /o^ cxp [Jp^ r(j/)(iy] dx. In general, the derivation of 
a stationary distribution of strategies under the effect of mu- 
tations cannot be solved analytically. Earlier approaches have 
approximated such distributions either by assuming weak mu- 
tations |33!3, where only transition probabilities between ab- 
sorbing states have to be considered, or by neglecting all mu- 
tations except those in the absorbing states I34ll . The distribu- 
tion given by Eq. ( I15> for large N is valid for arbitrary selec- 
tion intensity w, mutation rate u, and any payoff matrix. 



A. Neutral dynamics 

In the frequency dependent Moran process, often weak se- 
lection, w <C 1, is considered |13, 25]. For weak selection, 
the interactions are only corrections to a background of neu- 
tral dynamics lEl lssll . Therefore, we first discuss the station- 
ary distribution for w = 0, where the game has no influ- 
ence on the fitness. In this case, r(a;) is a gradient, which 
is equivalent to the assumption that the probability current 
vanishes everywhere, and the stationary distribution can be 
computed directly from Eq. ( II H . For d = 2, we find for 
the drift term a{x) = u{l — 2x) and for the diffusion term 
b{x) = {u{2x - 1)^ + 2x{l - x))/N. Both results are valid 
for the case in which the mutation rates from Ato B and vice 
versa are u. Hence, r(a;) is given by 
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, (m(2 + 7V) - l)(2a;- 1) 
'u{2x-l)'^ + 2x{l-x)' 



(16) 



r(a;) vanishes for the critical mutation rate Uc = 1/(2 + N). 
For u < lie, r(a;) has a minimum at x — 1/2. Accord- 
ingly, the stationary probability distribution has maxima at 
the boundaries when only few mutations per generation oc- 
cur On the other hand, for u > Uc there is a maximum at 
X = 1/2. Mutations lead the system away from the states 
X = and a; = 1 and the distribution is centered around the 
middle where both strategies have equal frequencies, see Fig. 
n For u = Uc we have r(a;) — and hence Eq. ( I15t predicts 





x=i/N 

FIG. 1 : For neutral selection, w = Q, three different scenarios occur 
for the stationary distribution of the system: For u < Uc = 1 / {2 + 
N), the system spends most of the time near the absorbing states, 
leading to maxima of the stationary distribution at i = and i = N, 
as shown for u = 0.005. For u = Uc ~ 0.01, the probability to 
leave the absorbing states due to mutations and the probability to 
reach them by random drift are approximately equal, leading to a 
uniform distribution. For u > Uc, there is on average more than one 
mutation per generation and the stationary distribution is centered 
around i = N/2. In all three cases, numerical simulations of the 
stationary distribution depicted by symbols agree very well with the 
theoretical result Eq. <I5> shown as lines (simulations are averages 
over 10** time steps, iV = 100). 

a uniform distribution. However, simulations show a distri- 
bution that decays slightly toward the boundaries for small N 
(cf. Fig. 0. Strictly speaking, a uniform distribution is ob- 
served only in the limit N ^ oo. 

Critical mutation rates can also be computed for general d 
in the case of symmetric mutations in which the mutation rate 
between two different states is u, i.e. qij = u+5ij{l — du). All 
elements of r(a;) vanish for u ^ Uc = + N), because 
2aj{x) ~ X]fe=i Hence, a uniform distribution 

is expected for li = in the limit of ^ oo. For smaller 
mutation rate u < u^, the distribution has maxima at the cor- 
ners of the simplex where only one strategy is present. For 
larger u > Uc, the distribution has one maximum in the inte- 
rior of the simplex. Fig.|2lillustrates these two different cases 
for neutral selection with three different types, d = i. 

B. Frequency dependent selection 

For w > 0, interactions lead to frequency dependent tran- 
sition probabilities and to corrections of the stationary distri- 
bution. Depending on the nature of the game and the strength 
of selection, significant differences from neutral evolution are 
possible. 

As examples, we consider two generic cases of 2 x 2 games, 
the Prisoner's Dilemma and the Snowdrift game. In the Pris- 
oner's Dilemma ll36ll37ll . two players choose simultaneously 
whether to cooperate or to defect. The cost of cooperation is 
c, while the benefit from cooperation is h. The interaction of 



FIG. 2: (Color online) For d = 3 strategies and neutral selection 
{w — 0), the stationary distribution in the strategy space spanned 
by the simplex S3 is shown encoded by a color scale, where bright 
colors indicate high values, (a) For mutation rates higher than the 
critical mutation rate Uc — 1/(3 + N), mutations drive the system 
away from the absorbing states at the corners of the simplex. Due to 
the symmetry of the system, the stationary distribution has a maxi- 
mum centered in the middle of the simplex where x = y — z. (b) 
If the mutation rates are smaller than than the critical mutation rate 
Uc, the system spends considerable time in the absorbing states at 
the corners. Occasionally, the edges are reached by mutations. How- 
ever, since the system typically reaches the corners again before the 
next mutation occurs, the stationary probability density in the inte- 
rior is very small. In both cases, the simulations of the stationary 
probability distribution shown here (averages over 10* time steps) 
do not deviate significantly from the analytical result Eq. <1 H : for 
u — 0.05 the maximal deviation is ~ 2% and for u = 0.005 it is 
~ 6% (population size — 60). 

a cooperator with a defector leads to a payoff of — c for the 
cooperator and b for the defector. When two cooperators in- 
teract, they both get the payoff b ~ c, whereas the interaction 
of two defectors leads to a payoff of 0. Thus, irrespective of 
the other player's move, it is always better to defect. How- 
ever, the reward for mutual cooperation would be the mutu- 
ally preferred outcome and hence the dilemma. The game is 
described by the payoff matrix 

^^/mn m,,\/b-c -c\ 
\ m2i m22 J \ b y ' 

The fitness is now a linear combination of the background 
fitness and the average payoff from interactions. For coopera- 
tors, the fitness becomes ttc = l~w+'w{(b—c){i—l)—c{N— 
i))/N , whereas defectors have fitness ttc = 1 — w + wbi/N. 
Since the fitness of defectors is always higher than the fitness 
of cooperators, ttd > ttc, any individual is better off not co- 
operating, despite the fact that mutual cooperation leads to 
a higher average payoff. Without mutations the majority of 
stochastic trajectories would therefore end up in the absorb- 
ing state with 100% defectors and only few trajectories end 
up in pure cooperation. For u < Uc = + 2), the station- 
ary distribution keeps a maximum at x = 0, but also a second 
(smaller) maximum at the other absorbing state x = 1, see 
Fig. |3l For u > Uc, mutations drive the system away from 
these states and a maximum appears for intermediate values 
of X. The position of this maximum is determined by the in- 
tensity of selection w and the payoff matrix. For increasing , 
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FIG. 3: Stationary distribution of strategies in a Prisoner's Dilemma 
under the influence of mutations and selection in a finite population. 
The distribution undergoes a qualitative change with increasing pop- 
ulation size. Main figure: For N — 50, the stationary distribution 
for w = 0.2 has a maximum at the Nash equilibrium a; = 0, where 
only defectors are present, as u is below the critical mutation rate 
Uc = 0.02. When selection is weak, it resembles the result from 
neutral selection (w — 0.01). Inset: When the population size is in- 
creased to A'^ = 10000, the mutation rate is above the critical value 
Uc ~ 10^* and the situation resembles the result expected from the 
deterministic description of the replicator-mutator equation: The sta- 
tionary distribution has a maximum at x* « 0.14 for w = 0.2, 
which is the stable fixed point of the replicator-mutator equation. 
For weak selection, w = 0.01, the stable fixed point is located at 
X* « 0.47 and the distribution around this fixed point is wider due 
to stronger fluctuations. In all cases, numerical simulations of the 
frequency dependent Moran process (symbols) agree well with the 
analytical solution Eq. <I5> for the stationary distribution depicted 
by lines {b = 1, c — 0.25, mutation rate u — 0.01, average over 10* 
time steps). 



the situation may change from u < Uc to u > Uc- Therefore, 
the stochastic effects arising from the finite population size 
are clearly visible for small N, whereas large N resembles a 
stationary distribution that is similar to the one expected from 
the deterministic replicator-mutator equation, see Sec.lVl 

As a second example, we consider the Snowdrift game llssll 
also known as "Hawk-Dove game" or "Chicken". In this 
game, two players again choose between cooperation and de- 
fection. While the payoffs for defectors are as in the Prisoner's 
Dilemma, cooperation becomes more favorable: The payoff 
of a cooperator against a defector is now b — c and therefore 
higher than the payoff of a defector against a defector. On the 
other hand, it is still lower than the payoff for mutual cooper- 
ation, b — c/2. The payoff matrix is given by 
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(18) 



Again, the payoff is a linear combination of the background 
fitness and the average payoff from interactions. In the Snow- 
drift game, rare strategies are always favored, resulting in a 
drift away from the absorbing states a; = and a; = 1. In- 
stead, the distribution is centered around a stable equilibrium 
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FIG. 4: The Snowdrift game exhibits a stable interior equilibrium 
and typically coexistence of the two different types is observed. Only 
for small mutation rates and weak selection, the absorbing states are 
reached. Therefore, the stationary distribution has a maximum in the 
interior For weak selection [w — 0.01) this maximum is close to 
X = 0.5. With stronger selection, the maximum moves toward the 
mixed equilibrium at x* « 0.856. Numerical simulations (symbols) 
are in good agreement with the theoretical result depicted by full 
lines (b = 1, c — 0.25, mutation rate u — 0.01, averages over 10* 
time steps. A*' = 200). 



in the interior. With decreasing intensity of selection, the dis- 
tribution around this equilibrium becomes wider. In the limit 
N oo and without mutations {u = 0), the distribution nar- 
rows and converges to x* ~ 2(6 — c)/{2b — c), the only sta- 
ble fixed point of the replicator dynamics. Finite size fluc- 
tuations lead to a distribution around this deterministic fixed 
point while mutations move the maximum of this distribution 
toward x = 0.5. As for the Prisoner's Dilemma, simulations 
of the stationary distribution agree well with the analytical re- 
sult Eq. ([B}, cf. Fig.H 



V. REPLICATOR-MUTATOR EQUATION 

In the limit of infinite population size, our framework es- 
tablishes a natural connection to the traditional determinis- 
tic description of coevolutionary dynamics. The determin- 
istic nature of this limit arises from the fact that the diffu- 
sion matrix bkj{x) and thus the finite size fluctuations vanish 
with I/^/N as N oo. Hence, the Langevin equation be- 
comes deterministic and the dynamics is given by the drift 
term, x^ = J2k=i i'^jk ~ Tkj). Since the mutation rate per 
individual is fixed, the effect of mutations prevails in the limit 
of iV ^ oo. In this case, Eq. Q incorporating mutation and 
selection can be written as 



Xk 
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Xk- 



(19) 



This is the replicator-mutator equation corresponding to the 
adjusted replicator dynamics 0]. So far, replicator-mutator 
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equations iHH [l2il have not been derived from a microscopic 
birth-death process. They are the natural extension of the de- 
terministic replicator equation fl(ill that take into account 
mutations between strategies. Without mutations, qjk = Sjk, 
types that have a fitness above the average fitness in the pop- 
ulation increase in abundance, while types with a fitness be- 
low average decrease in abundance. Mutations interfere with 
this selection dynamics: Only types that have a fitness above 
the average and that have a sufficiently small mutation rate 
increase in frequency. Hence, the stable fixed points of the 
system represent a balance between selection dynamics and 
mutation dynamics. 



VI. LOCAL UPDATE PROCESS 

So far, we have only considered the frequency dependent 
Moran process. However, our framework can also be applied 
to other coevolutionary processes. As an example, we gener- 
alize the local update process f29'| to an arbitrary number of 
strategies. In the local update process, an individual 1 is cho- 
sen at random and compares its payoff to another randomly 
chosen individual 2. With probability 
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2" Att 



(20) 



1 adapts the strategy of 2 (where the constant Att is the max- 
imal possible payoff difference and w denotes the intensity of 
selection) |39]. For d different types, the transition probabili- 
ties are given by 



(21) 



From the transition probabilities, we can calculate the drift 
term ak{x) = (w/Att) Xk {TTk{x) — </>), which is indepen- 
dent of N. The diffusion term is independent of the payoffs 
and simplifies to bjk{x) = Xj {Sjk — Xk) /N. Thus, up to a 
constant rescaling of time, we recover the standard replicator 
dynamics for d types in the limit N ^ oo. 



Xk = Xk (TTkix) - (/)) . 



(22) 



Up to a dynamical rescaling of time, Eq. i22i is identical to 
Eq. ( I19> without mutations. The natural extension of Eq. (I22> 
that includes mutations (while retaining the time scale) is the 
standard replicator-mutator equation 



ing tt 
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XjTTj{x)qjk-Xk(l> 
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Xj TTj (x) {xiQjk - XkQkl ) • 

(23) 

For general qjk, it is not possible to define a microscopic pro- 
cess in which the transition probabilities only depend on pay- 
off differences that leads to this differential equation. Only 
for the special case of vanishing mutations qjk — Sjk, this is 
possible. 

A microscopic mechanism involving spontaneous muta- 
tions has been proposed in ll40ll . This yields the transition 



rates 141L 

, X /I W TTiix) - 'Kk(x)\ 

Tkj (x) = XkX, + + a^fcgfcj . (24) 

As shown in j40ll . the limit iV ^ oo results in the differential 
equation 



Xk 



Att 



Xk {T^kix) - 0) + ^ {xjQjk - XkQkj) ■ (25) 



This selection mutation equation reduces to the standard repli- 
cator dynamics for qjk = Sjk Em. It describes selection 
under the influence of spontanteous mutations. Such sponta- 
neous mutations can also be incorporated into other coevolu- 
tionary processes, leading to the same additive mutation term 

ELi {Xj(ljk - XkQkj)- 



VIL DISCUSSION 

Traditional descriptions of co-evolutionary systems rely on 
the assumption that populations are sufficiently large such that 
they can safely be considered infinite. In this case, accidental 
fixation of types with low fitness is not possible. For the more 
realistic case of finite populations, only few analytical meth- 
ods exist that are often restricted to two types. Hence, these 
systems often rely on simulations which are very illuminative 
for special cases, but give little insight into the general dynam- 
ics of these systems. Here, we have shown how to extend the 
replicator dynamics in order to account for fluctuations arising 
from finite population sizes for any number of strategies and 
arbitrary mutation rates. When allowing for mutations, homo- 
geneous states in which only one type is present are no longer 
absorbing boundaries and, instead, the system converges to a 
stationary distribution. This stationary distribution is derived 
using standard methods from statistical physics and depends 
on the population size, the mutation rate and the interaction 
parameters, i.e., the game. The fluctuation term allows for the 
calculation of corrections to the replicator dynamics for finite 
populations and is an important step in comparing individ- 
ual based coevolutionary systems with deterministic dynam- 
ics. Mutations are not a source of fluctuations, but provide a 
mechanism that leads to a mixing of different types. However, 
fluctuations do not arise solely from finite populations. In- 
stead, there are additional sources of noise that might have to 
be incorporated in other ways, such as noise from heterogene- 
ity of agents or external noise that can influence interactions. 

Our approach can be applied to any number of strategies d, 
but is based on the assumption that the population size is suf- 
ficiently large such that population densities can be approxi- 
mated by smooth functions. In high dimensional spaces, such 
a continuum description may not always be appropriate 1 43 ] . 

Starting from microscopic descriptions for individual based 
coevolutionary processes with an arbitrary number of strate- 
gies, we have derived a replicator mutator equation with ad- 
ditional first order corrections that describe stochastic effects 
arising from finite populations. 
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APPENDIX A: DERIVATION OF THE FOKKER-PLANCK 
EQUATION FOR AN ARBITRARY NUMBER OF 
STRATEGIC TYPES 

For d strategies, the Master equation describing the change 
of the probability P^(ii, .., id) that the system is in focal state 
{ii, ..,id) is given by 



P 



t+At 



{ii,..,id) = P^{ii,-;id) (Al) 

d 

3,k=l 

X Tkj{ii, ...,ik + l, ■■■,id) 

d 

-At^ P^{ii,...,id)Tkj{ii,...,id). 



Other states. The probability Tkj{ii, ■■■,id) that one individ- 
ual of type k is replaced by another individual of type j does 
not have to be specified. Special cases are given by Eqs. Q 
and i24l . Next, we introduce the notation xi = ii/N and 
t = t/N and the probability density p{x; t) ~ NP*'{x). For 
simplicity, we use the abbreviation (7^, fc+) = (xi, ...Xj — 
1/iV, + 1/iV, ...Xd) for states neighboring the focal state 
X ~ (ii, . . . , id)/N and obtain 



p{x]t + At) - p{x;t) 
NAt 



pu-,k+;t)n^{j-,k+) 

],k=i 
d 

- p{x;t)Tk,{x), (A2) 

3,k=l 



The first term is the probability that during a short time inter- 
val At the state remains unchanged, the second term is the 
influx from other states and the third term is the outflux to 



For N ^ 1, the left hand side can be approximated by 
p{x; t) /N . In order to expand the right hand side of Eq. JA2t 
for N ':$> 1, we rewrite the first term as 



J 



Y pir,k+;t)Tk,{r,k+)^ Y 
j.k=i j,k=i 



p{j-,k+;t)n,U-,k+) + p{j+, fc-; t)T,fc(j+, k-) 



The expansion of the probability density at a state {j^ ,k^) at time t can be written as 



p{j^,k^;t)^p{x;t)±- 



dp{x-t) dp{x;t) 
dxi 



dxk 



and similarly for the transition probabilities 



1 



Tkji3^,k^)^Tk,ix)±- 



dTkj{x) dTkjix) 



dxi 



dxk 



2N^ 



27V2 



d'p{x;t) ^d^p{x-t) , d'p{x-t) 



dx] 



dxjdxk 



dxl 



d^Tk,{x) ^d^Tk,{x) , d^Tk,{x) 



dx^ 



dxjdxk 



dxl 



(A3) 



(A4) 



(A5) 



Let us consider the terms on the right hand side of Eq. iA2\ depending on their order in 1 /N. Obviously, the terms that are 
independent of N cancel. Up to we find 



1 

27V 



E 



pix;t) ( ^ - 'JkM.) + (9_Me1 _ 9J^) ^TkAx) C'^^-^'^ '^^^-''^ 



+Tjk{x) 



dxk dxj 
f dp{x;t) dp{x;t) 



dxj 



V dxk 



dxi 



V dxj 



dxk 



(A6) 



k=l 



i=i 



Since 1 /N appearing on both sides of Eq. iA2\ cancels, this determines the drift term that is independent of N. All other terms 
vanish for N 00. The first correction for finite population sizes is the diffusion term resulting from terms in the expansion of 
the right hand side of Eq. (IA2> that are of order 1/N'^. For these terms, we find 



1 

m E 



2Ar2 



j,k=i L 



Tkj{x) + T,k{x) / 92 



dx] 



92 92 \ / Q2 Q2 

2 „ „ + — 1 pix; t) + p{x- 1)\ —-2- 



92 \ Tkj{x)+T,k{x) 



dxidxk 



dxl 



dx] dxjdxk dx\ 



dp{x-t) dTkj{x) dp{x-t)dTkj{x) dp{x-t)dTk,{x) dp{x-t)dTk,{x) dp{x;t) dT.kix) dp{x-t)dTjk{x) 



dxk 



dxk 



dxk 



dxi 



dp{x-t)dTjk{x) , dp{x-t)dTjk{s 



dx-i 



dxk 



dxj 



dxi 



dxj 



2iV2 



dxi 



dxi 



E 

fe,i=i 



dx^. dxjdxk 



dxj dxk dxk 

p{x-t) {Tk,[x)+T,k{x)) 



dxk 



dxj 
(A7) 
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So far, we have not used the normaHzation condition, and a diffusion matrix bji defined as 



= 1. Since is fully determined by xj with 



I < j < d, we eliminate the dependence on Xd in the transi- 
tion probabilities Tjk and in the probability density p. Hence, 
all derivatives with respect to Xd vanish. Applying this and 
combining (IA6> and iA7\ with the expansion of the left hand 
side of Eq. (IA2> . we obtain the Fokker-Planck equation 



fc=l 

with drift vector 



2 dxkdxj 

Jtf^ — -I- 



p{x)bjk{x) 
(A8) 



(A9) 



hjk{x) 



N 



-T,k{x)~ [x) + 5,kY\ Tji (x) + Ti, {x) 



1=1 



(AlO) 

which is symmetric, bjk{x) = bkj{x). Note that the drift 
vector of a system with d strategies has — 1 components, 
/c = 1 , . . . , (i — 1 and that the diffusion matrix is of dimension 
1. 
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